library(qualtRics)
library(stargazer)
library(ggplot2)
library(dplyr)
library(httr)
library(data.table)
library(vtable)
library(stargazer)

setwd(dirname(rstudioapi::getSourceEditorContext()$path))

###Importing and Processing Survey Data###
# EXPERIMENTS 1 and 2 – We ran these on the same day and manipulate the variables using the same code.

# Import Survey data from Experiments 1 and 2 and correct variable names
df1 <- readRDS("DataFiles\\Experiment1_Iran_Raw.rds")
df2 <- readRDS("DataFiles\\Experiment2_Qatar_Raw.rds")
df3 <- rbindlist(list(df1, df2), fill = TRUE)
names(df3) <- gsub(":", "", names(df3))
names(df3) <- make.names(names(df3), unique=TRUE)
###Attention Checks###

# Drop based on first attention check
df <- df3[
  !is.na(df3$AC1_1) & 
    !is.na(df3$AC1_2) & 
    df3$Status != "Survey Preview" & df3$Finished == "TRUE",]

# Convert string NAs to Actual NAs
df[] <- lapply(df, function(x) {
  is.na(levels(x)) <- levels(x) == "NA"
  x
})

# Drop based on second attention check
df = df[df$i.AC2 == "The United States and Iran" | df$q.AC2 == "The United States and Qatar",]

# Keep only respondents who finished and were from Lucid (no previews)
df <- df[df$Status == "IP Address",]
df <- df[!is.na(df$Status),]

# Check what percentage of respondents got the manipulation check correct
df$AC3 = ifelse(  df$i.AC3 == "Iran claimed responsibility for the attack." & !is.na(df$i.overt) |
                  df$i.AC3 == "Iran denied involvement in the attack." & !is.na(df$i.denial) |
                  df$q.AC3 == "Qatar claimed responsibility for the attack." & !is.na(df$q.overt) |
                  df$q.AC3 == "Qatar denied involvement in the attack." & !is.na(df$q.denial),1,0)
perc3 = sum(df$AC3)/length(df$AC3)
  

# Create denial and adversary variables
df$denial = ifelse(!is.na(df$q.denial) | !is.na(df$i.denial), 1,0)
df$adversary = ifelse(!is.na(df$i.AC2), 1,0)

###Combining Variables###

df = df %>% mutate(escalation_1 = coalesce(q.escalation_1, i.escalation_1))
df = df %>% mutate(escalation_2 = coalesce(q.escalation_2, i.escalation_2))
df = df %>% mutate(escalation_3 = coalesce(q.escalation_3, i.escalation_3))
df = df %>% mutate(escalation_4 = coalesce(q.escalation_4, i.escalation_4))

df = df %>% mutate(reputation_1 = coalesce(q.reputation_1, i.reputation_1))
df = df %>% mutate(reputation_2 = coalesce(q.reputation_2, i.reputation_2))
df = df %>% mutate(reputation_3 = coalesce(q.reputation_3, i.reputation_3))

df = df %>% mutate(insulting = coalesce(q.emotions, i.emotions))

df = df %>% mutate(ambiguity = coalesce(q.ambiguity, i.ambiguity))

###Creating Indices for Variables###

# Factor escalation variables
df$escalation_1 = factor(df$escalation_1, levels = c(
  "Oppose strongly",
  "Oppose somewhat",
  "Neither favor nor oppose",
  "Favor somewhat",
  "Favor strongly"
))

df$escalation_2 = factor(df$escalation_2, levels = c(
  "Oppose strongly",
  "Oppose somewhat",
  "Neither favor nor oppose",
  "Favor somewhat",
  "Favor strongly"
))

df$escalation_3 = factor(df$escalation_3, levels = c(
  "Oppose strongly",
  "Oppose somewhat",
  "Neither favor nor oppose",
  "Favor somewhat",
  "Favor strongly"
))

df$escalation_4 = factor(df$escalation_4, levels = c(
  "Oppose strongly",
  "Oppose somewhat",
  "Neither favor nor oppose",
  "Favor somewhat",
  "Favor strongly"
))

# Scale and combine escalation variables
df$esca_scaled = (scale(as.numeric(df$escalation_1)) + 
                    scale(as.numeric(df$escalation_2)) + 
                    scale(as.numeric(df$escalation_3)) + 
                    scale(as.numeric(df$escalation_4)))/4
df$esca_scaled = (df$esca_scaled - min(df$esca_scaled, na.rm = TRUE))
df$esca_scaled = df$esca_scaled/max(df$esca_scaled,na.rm = TRUE)*100

# Weighted average version
df$esca_scaled_wa = (scale(as.numeric(df$escalation_1))*1 + 
                    scale(as.numeric(df$escalation_2))*2 + 
                    scale(as.numeric(df$escalation_3))*3 + 
                    scale(as.numeric(df$escalation_4)))*4/(4+3+2+1)
df$esca_scaled_wa = (df$esca_scaled_wa - min(df$esca_scaled_wa, na.rm = TRUE))
df$esca_scaled_wa = df$esca_scaled_wa/max(df$esca_scaled_wa,na.rm = TRUE)*100

# Create individual escalation variables for each answer
df$war = df$escalation_4
df$airstrike = df$escalation_3
df$sanctions = df$escalation_2
df$diplomacy = df$escalation_1

# Scale Military Assertiveness
df$MA = (6 - as.numeric(df$MA1)) + as.numeric(df$MA2) + (6 - as.numeric(df$MA3))
df$MA_scaled = (scale(6 - as.numeric(df$MA1))) + 
  scale(as.numeric(df$MA2)) + 
  scale((6 - as.numeric(df$MA3)))/3

# Scale National Chauvinism
df$NC = (as.numeric(df$NC1) + as.numeric(df$NC2))
df$NC_scaled = scale(as.numeric(df$NC1) + scale(as.numeric(df$NC2)))/2

# Factor reputation variables
df$reputation_1 = factor(df$reputation_1, levels = c(
  "Almost no chance",
  "25% change",
  "50-50 chance",
  "75% chance",
  "Nearly 100% certain"
))

df$reputation_2 = factor(df$reputation_2, levels = c(
  "Almost no chance",
  "25% change",
  "50-50 chance",
  "75% chance",
  "Nearly 100% certain"
))

df$reputation_3 = factor(df$reputation_3, levels = c(
  "Almost no chance",
  "25% change",
  "50-50 chance",
  "75% chance",
  "Nearly 100% certain"
))

# Scale and combine reputation variables
df$reputation_scaled = (scale(as.numeric(df$reputation_1)) + scale(as.numeric(df$reputation_2)) + scale(as.numeric(df$reputation_3)))/3
df$reputation_scaled = (df$reputation_scaled - min(df$reputation_scaled, na.rm = TRUE))
df$reputation_scaled = df$reputation_scaled/max(df$reputation_scaled,na.rm = TRUE)*100

# Convert dispositional controls to numeric values
df$GovTrust = 6-as.numeric(df$GovTrust)
df$NewsTrust = 6-as.numeric(df$NewsTrust)
df$IntTrust = 6-as.numeric(df$IntTrust)
df$Military.Service = 2-as.numeric(df$Military.Service)
df$Military.Service[df$Military.Service ==-1] = NA
df$Read.FP = 6-as.numeric(df$Read.FP)

###Cleaning demographic controls###

# Make male variable 0 and 1
df$male = df$gender
df$male[df$male == 2] = 0

# Change missing hhi data to NA
df$hhi[df$hhi == -3105] = NA

# Create hispanic indicator variable
df$hispanic[df$hispanic == 15] = NA
df$hispanic = ifelse(df$hispanic > 1, 1, 0)

# Create white and black indicator variables
df$white = ifelse(df$ethnicity == 1 & df$hispanic == 0, 1, 0)
df$black = ifelse(df$ethnicity == 2, 1, 0)

# Change missing education data to 0
df$education[df$education == -3105] = 0

# Create republican and democrat indicator variables (independent comparison)
df$republican = ifelse(
  df$political_party == 9 | 
    df$political_party == 10 | 
    df$political_party == 5 | 
    df$political_party == 8, 1, 0)
df$democrat = ifelse(
  df$political_party == 1 | 
    df$political_party == 2 | 
    df$political_party == 3 | 
    df$political_party == 6, 1, 0)

###Final Dataset Export###

colnames(df)[colnames(df) == 'ambiguity'] <- 'certainty'
df_res = df %>% dplyr::select(denial, adversary, esca_scaled, esca_scaled_wa, MA_scaled, GovTrust, 
                              NewsTrust, IntTrust, NC_scaled, Military.Service, Read.FP,
                              reputation_scaled, certainty, insulting, war, airstrike, sanctions, 
                              diplomacy, age, male, hhi, white, black, education, republican, democrat, AC3)
df_res[] <- lapply(df_res, as.numeric)
df_res$ResponseId = df$ResponseId

# Experiment 1 (Iran) Data
df_i = df_res[df_res$adversary==1]
df_i$adversary = NULL
write.csv(df_i, "DataFiles\\Experiment1_Iran.csv", row.names = FALSE)

# Experiment 2 (Qatar) Data
df_q = df_res[df_res$adversary==0]
df_q$adversary = NULL
write.csv(df_q, "DataFiles\\Experiment2_Qatar.csv", row.names = FALSE)

###Bar Plot Iran####

# Subset for Iran

# Calculate means of relevant variables
means = c(mean(df_i$war[df_i$denial == 1],na.rm = TRUE), 
          mean(df_i$airstrike[df_i$denial == 1],na.rm = TRUE),
          mean(df_i$sanctions[df_i$denial == 1],na.rm = TRUE), 
          mean(df_i$diplomacy[df_i$denial == 1],na.rm = TRUE),
          mean(df_i$war[df_i$denial == 0],na.rm = TRUE), 
          mean(df_i$airstrike[df_i$denial == 0],na.rm = TRUE),
          mean(df_i$sanctions[df_i$denial == 0],na.rm = TRUE), 
          mean(df_i$diplomacy[df_i$denial == 0],na.rm = TRUE)
          )

# Calculate standard deviations of relevant variables
sds = c(sd(df_i$war[df_i$denial == 1],na.rm = TRUE), 
          sd(df_i$airstrike[df_i$denial == 1],na.rm = TRUE),
          sd(df_i$sanctions[df_i$denial == 1],na.rm = TRUE), 
          sd(df_i$diplomacy[df_i$denial == 1],na.rm = TRUE),
          sd(df_i$war[df_i$denial == 0],na.rm = TRUE), 
          sd(df_i$airstrike[df_i$denial == 0],na.rm = TRUE),
          sd(df_i$sanctions[df_i$denial == 0],na.rm = TRUE), 
          sd(df_i$diplomacy[df_i$denial == 0],na.rm = TRUE)
)

# Calculate N of relevant variables
ns = c(length(df_i$war[df_i$denial == 1]), 
        length(df_i$airstrike[df_i$denial == 1]),
        length(df_i$sanctions[df_i$denial == 1]), 
        length(df_i$diplomacy[df_i$denial == 1]),
        length(df_i$war[df_i$denial == 0]), 
        length(df_i$airstrike[df_i$denial == 0]),
        length(df_i$sanctions[df_i$denial == 0]), 
        length(df_i$diplomacy[df_i$denial == 0])
)

# Create data frame with labels and data
treat = c("Covert", "Covert","Covert","Covert", "Overt", "Overt", "Overt", "Overt")
type = c("War", "Airstrike", "Sanctions", "Condemn", "War", "Airstrike", "Sanctions", "Condemn")
type = factor(type, levels = c(
  "Condemn",
  "Sanctions",
  "Airstrike",
  "War"
))
forgraph = data.frame(means, treat, type, sds, ns)
forgraph$error = qt(0.975, df=ns-1)*sds/sqrt(ns)

# Figure 2: Support for Escalation Options by Treatment Condition in Experiment 1
ggplot(aes(x = type, y = means, fill = treat), data = forgraph, group = factor(type)) +
  geom_bar(stat = "identity", 
           position = position_dodge(width = 0.9)) +
  geom_errorbar(aes(x=type, ymin=means-error, ymax=means+error),
                width = 0.25,
                position = position_dodge(width = 0.9)) +
  xlab("") +
  ylab("") +
  scale_y_discrete(limits = c(1,2,3,4,5),
                   labels = c("Strongly\nOppose (1)", 
                              "Somewhat\nOppose (2)",
                              "Neutral (3)", 
                              "Somewhat\nFavor (4)", 
                              "Strongly\nFavor (5)"),
                   expand = expansion(add = c(0,1))) +
  scale_fill_discrete(name = "", labels = c("Denial", "Overt")) +
  theme(axis.text.y = element_text(angle = 0, size = 12.5),
        axis.text.x = element_text(size = 12.5),
        legend.text = element_text(size = 12.5)) 
ggsave("Tables&Figures\\avg_escalation_iran.png", width = 6, height = 4, unit = "in")

###Bar Plot Qatar###

# Subset for Qatar

# Calculate means of relevant variables
means = c(mean(df_q$war[df_q$denial == 1],na.rm = TRUE), 
          mean(df_q$airstrike[df_q$denial == 1],na.rm = TRUE),
          mean(df_q$sanctions[df_q$denial == 1],na.rm = TRUE), 
          mean(df_q$diplomacy[df_q$denial == 1],na.rm = TRUE),
          mean(df_q$war[df_q$denial == 0],na.rm = TRUE), 
          mean(df_q$airstrike[df_q$denial == 0],na.rm = TRUE),
          mean(df_q$sanctions[df_q$denial == 0],na.rm = TRUE), 
          mean(df_q$diplomacy[df_q$denial == 0],na.rm = TRUE)
)

# Calculate standard deviations of relevant variables
sds = c(sd(df_q$war[df_q$denial == 1],na.rm = TRUE), 
        sd(df_q$airstrike[df_q$denial == 1],na.rm = TRUE),
        sd(df_q$sanctions[df_q$denial == 1],na.rm = TRUE), 
        sd(df_q$diplomacy[df_q$denial == 1],na.rm = TRUE),
        sd(df_q$war[df_q$denial == 0],na.rm = TRUE), 
        sd(df_q$airstrike[df_q$denial == 0],na.rm = TRUE),
        sd(df_q$sanctions[df_q$denial == 0],na.rm = TRUE), 
        sd(df_q$diplomacy[df_q$denial == 0],na.rm = TRUE)
)

# Calculate N of relevant variables
ns = c(length(df_q$war[df_q$denial == 1]), 
       length(df_q$airstrike[df_q$denial == 1]),
       length(df_q$sanctions[df_q$denial == 1]), 
       length(df_q$diplomacy[df_q$denial == 1]),
       length(df_q$war[df_q$denial == 0]), 
       length(df_q$airstrike[df_q$denial == 0]),
       length(df_q$sanctions[df_q$denial == 0]), 
       length(df_q$diplomacy[df_q$denial == 0])
)

# Create data frame with labels and data
treat = c("Covert", "Covert","Covert","Covert", "Overt", "Overt", "Overt", "Overt")
type = c("War", "Airstrike", "Sanctions", "Condemn", "War", "Airstrike", "Sanctions", "Condemn")
type = factor(type, levels = c(
  "Condemn",
  "Sanctions",
  "Airstrike",
  "War"
))
forgraph = data.frame(means, treat, type, sds, ns)
forgraph$error = qt(0.975, df=ns-1)*sds/sqrt(ns)

# Figure 4: Support for Escalation Options by Treatment Condition in Experiment 2
ggplot(aes(x = type, y = means, fill = treat), data = forgraph, group = factor(type)) +
  geom_bar(stat = "identity", 
           position = position_dodge(width = 0.9)) +
  geom_errorbar(aes(x=type, ymin=means-error, ymax=means+error),
                width = 0.25,
                position = position_dodge(width = 0.9)) +
  xlab("") +
  ylab("") +
  ggtitle("Qatar: Average Preference per Response Option") +
  scale_y_discrete(limits = c(1,2,3,4,5),
                   labels = c("Strongly\nOppose (1)", 
                              "Somewhat\nOppose (2)", 
                              "Neutral (3)", 
                              "Somewhat\nFavor (4)", 
                              "Strongly\nFavor (5)"),
                   expand = expansion(add = c(0,1))) +
  scale_fill_discrete(name = "", labels = c("Denial", "Overt")) +
  theme(axis.text.y = element_text(angle = 0, size = 12.5),
        axis.text.x = element_text(size = 12.5),
        legend.text = element_text(size = 12.5)) 
ggsave("Tables&Figures\\avg_escalation_qatar.png", width = 6, height = 4, unit = "in")


# Experiment 3 (Iran, More Certainty) --------------------------------------
# We ran this experiment later and cleaned the data separately.

# Import Survey and correct variable names
df <- readRDS("DataFiles\\Experiment3_Iran_Raw.rds")
names(df) <- gsub(":", "", names(df))
names(df) <- make.names(names(df), unique=TRUE)
df = as.data.frame(df)
###Attention Checks###

# Drop based on first attention check
df <- df[
  !is.na(df$AC1_1) & 
    !is.na(df$AC1_2) & 
    df$Status != "Survey Preview" & df$Finished == "TRUE",]

# Convert string NAs to Actual NAs
df[] <- lapply(df, function(x) {
  is.na(levels(x)) <- levels(x) == "NA"
  x
})

# Drop based on second attention check
df = df[df$i.AC2 == "The United States and Iran",]

# Keep only respondents who finished and were from Lucid (no previews)
df <- df[df$Status == "IP Address",]
df <- df[!is.na(df$Status),]

# Check what percentage of respondents got the manipulation check correct
df$AC3 = ifelse(  df$i.AC3 == "Iran claimed responsibility for the attack." & !is.na(df$i.overt) |
                    df$i.AC3 == "Iran denied involvement in the attack." & !is.na(df$i.denial),1,0)
perc3 = sum(df$AC3)/length(df$AC3)

# Create denial and adversary variables
df$denial = ifelse(!is.na(df$i.denial), 1,0)
df$adversary = ifelse(!is.na(df$i.AC2), 1,0)

###Combining Variables###

df = df %>% mutate(escalation_1 = coalesce(i.escalation_1))
df = df %>% mutate(escalation_2 = coalesce(i.escalation_2))
df = df %>% mutate(escalation_3 = coalesce(i.escalation_3))
df = df %>% mutate(escalation_4 = coalesce(i.escalation_4))

df = df %>% mutate(reputation_1 = coalesce(i.reputation_1))
df = df %>% mutate(reputation_2 = coalesce(i.reputation_2))
df = df %>% mutate(reputation_3 = coalesce(i.reputation_3))

df = df %>% mutate(insulting = coalesce(i.emotions))

df = df %>% mutate(ambiguity = coalesce(i.ambiguity))

###Creating Indices for Variables###

# Factor escalation variables
df$escalation_1 = factor(df$escalation_1, levels = c(
  "Oppose strongly",
  "Oppose somewhat",
  "Neither favor nor oppose",
  "Favor somewhat",
  "Favor strongly"
))

df$escalation_2 = factor(df$escalation_2, levels = c(
  "Oppose strongly",
  "Oppose somewhat",
  "Neither favor nor oppose",
  "Favor somewhat",
  "Favor strongly"
))

df$escalation_3 = factor(df$escalation_3, levels = c(
  "Oppose strongly",
  "Oppose somewhat",
  "Neither favor nor oppose",
  "Favor somewhat",
  "Favor strongly"
))

df$escalation_4 = factor(df$escalation_4, levels = c(
  "Oppose strongly",
  "Oppose somewhat",
  "Neither favor nor oppose",
  "Favor somewhat",
  "Favor strongly"
))

# Scale and combine escalation variables
df$esca_scaled = (scale(as.numeric(df$escalation_1)) + 
                    scale(as.numeric(df$escalation_2)) + 
                    scale(as.numeric(df$escalation_3)) + 
                    scale(as.numeric(df$escalation_4)))/4
df$esca_scaled = (df$esca_scaled - min(df$esca_scaled, na.rm = TRUE))
df$esca_scaled = df$esca_scaled/max(df$esca_scaled,na.rm = TRUE)*100

# Weighted average version
df$esca_scaled_wa = (scale(as.numeric(df$escalation_1))*1 + 
                       scale(as.numeric(df$escalation_2))*2 + 
                       scale(as.numeric(df$escalation_3))*3 + 
                       scale(as.numeric(df$escalation_4)))*4/(4+3+2+1)
df$esca_scaled_wa = (df$esca_scaled_wa - min(df$esca_scaled_wa, na.rm = TRUE))
df$esca_scaled_wa = df$esca_scaled_wa/max(df$esca_scaled_wa,na.rm = TRUE)*100

# Create individual escalation variables for each answer
df$war = df$escalation_4
df$airstrike = df$escalation_3
df$sanctions = df$escalation_2
df$diplomacy = df$escalation_1

# Scale Military Assertiveness
df$MA = (6 - as.numeric(df$MA1)) + as.numeric(df$MA2) + (6 - as.numeric(df$MA3))
df$MA_scaled = (scale(6 - as.numeric(df$MA1))) + 
  scale(as.numeric(df$MA2)) + 
  scale((6 - as.numeric(df$MA3)))/3

# Scale National Chauvinism
df$NC = (as.numeric(df$NC1) + as.numeric(df$NC2))
df$NC_scaled = scale(as.numeric(df$NC1) + scale(as.numeric(df$NC2)))/2

# Factor reputation variables
df$reputation_1 = factor(df$reputation_1, levels = c(
  "Almost no chance",
  "25% change",
  "50-50 chance",
  "75% chance",
  "Nearly 100% certain"
))

df$reputation_2 = factor(df$reputation_2, levels = c(
  "Almost no chance",
  "25% change",
  "50-50 chance",
  "75% chance",
  "Nearly 100% certain"
))

df$reputation_3 = factor(df$reputation_3, levels = c(
  "Almost no chance",
  "25% change",
  "50-50 chance",
  "75% chance",
  "Nearly 100% certain"
))

# Scale and combine reputation variables
df$reputation_scaled = (scale(as.numeric(df$reputation_1)) + scale(as.numeric(df$reputation_2)) + scale(as.numeric(df$reputation_3)))/3
df$reputation_scaled = (df$reputation_scaled - min(df$reputation_scaled, na.rm = TRUE))
df$reputation_scaled = df$reputation_scaled/max(df$reputation_scaled,na.rm = TRUE)*100

# Convert dispositional controls to numeric values
df$GovTrust = 6-as.numeric(df$GovTrust)
df$NewsTrust = 6-as.numeric(df$NewsTrust)
df$IntTrust = 6-as.numeric(df$IntTrust)
df$Read.FP = 6-as.numeric(df$Read.FP)

###Cleaning demographic controls###

# Make male variable 0 and 1
df$male = df$gender
df$male[df$male == 2] = 0

# Change missing hhi data to NA
df$hhi[df$hhi == -3105] = NA

# Create hispanic indicator variable
df$hispanic[df$hispanic == 15] = NA
df$hispanic = ifelse(df$hispanic > 1, 1, 0)

# Create white and black indicator variables
df$white = ifelse(df$ethnicity == 1 & df$hispanic == 0, 1, 0)
df$black = ifelse(df$ethnicity == 2, 1, 0)

# Change missing education data to 0
df$education[df$education == -3105] = 0

# Create republican and democrat indicator variables (independent comparison)
df$republican = ifelse(
  df$political_party == 9 | 
    df$political_party == 10 | 
    df$political_party == 5 | 
    df$political_party == 8, 1, 0)
df$democrat = ifelse(
  df$political_party == 1 | 
    df$political_party == 2 | 
    df$political_party == 3 | 
    df$political_party == 6, 1, 0)

###Final Data Export Experiment 3###

colnames(df)[colnames(df) == 'ambiguity'] <- 'certainty'
df_res = df %>% dplyr::select(denial, esca_scaled, esca_scaled_wa, MA_scaled, GovTrust, 
                              NewsTrust, IntTrust, NC_scaled, Read.FP,
                              reputation_scaled, certainty, insulting, war, airstrike, sanctions, 
                              diplomacy, age, male, hhi, white, black, education, republican, democrat, AC3)
df_res[] <- lapply(df_res, as.numeric)
df_res$ResponseId = df$ResponseId
write.csv(df_res, "DataFiles\\Experiment3_Iran.csv", row.names = FALSE)

###Bar Plot Iran 3####

# Subset for Iran
# Calculate means of relevant variables
means = c(mean(df_res$war[df_res$denial == 1],na.rm = TRUE), 
          mean(df_res$airstrike[df_res$denial == 1],na.rm = TRUE),
          mean(df_res$sanctions[df_res$denial == 1],na.rm = TRUE), 
          mean(df_res$diplomacy[df_res$denial == 1],na.rm = TRUE),
          mean(df_res$war[df_res$denial == 0],na.rm = TRUE), 
          mean(df_res$airstrike[df_res$denial == 0],na.rm = TRUE),
          mean(df_res$sanctions[df_res$denial == 0],na.rm = TRUE), 
          mean(df_res$diplomacy[df_res$denial == 0],na.rm = TRUE)
)

# Calculate standard deviations of relevant variables
sds = c(sd(df_res$war[df_res$denial == 1],na.rm = TRUE), 
        sd(df_res$airstrike[df_res$denial == 1],na.rm = TRUE),
        sd(df_res$sanctions[df_res$denial == 1],na.rm = TRUE), 
        sd(df_res$diplomacy[df_res$denial == 1],na.rm = TRUE),
        sd(df_res$war[df_res$denial == 0],na.rm = TRUE), 
        sd(df_res$airstrike[df_res$denial == 0],na.rm = TRUE),
        sd(df_res$sanctions[df_res$denial == 0],na.rm = TRUE), 
        sd(df_res$diplomacy[df_res$denial == 0],na.rm = TRUE)
)

# Calculate N of relevant variables
ns = c(length(df_res$war[df_res$denial == 1]), 
       length(df_res$airstrike[df_res$denial == 1]),
       length(df_res$sanctions[df_res$denial == 1]), 
       length(df_res$diplomacy[df_res$denial == 1]),
       length(df_res$war[df_res$denial == 0]), 
       length(df_res$airstrike[df_res$denial == 0]),
       length(df_res$sanctions[df_res$denial == 0]), 
       length(df_res$diplomacy[df_res$denial == 0])
)

# Create data frame with labels and data
treat = c("Covert", "Covert","Covert","Covert", "Overt", "Overt", "Overt", "Overt")
type = c("War", "Airstrike", "Sanctions", "Condemn", "War", "Airstrike", "Sanctions", "Condemn")
type = factor(type, levels = c(
  "Condemn",
  "Sanctions",
  "Airstrike",
  "War"
))
forgraph = data.frame(means, treat, type, sds, ns)
forgraph$error = qt(0.975, df=ns-1)*sds/sqrt(ns)

# Figure 6: Support for Escalation Options by Treatment Condition in Experiment 3
ggplot(aes(x = type, y = means, fill = treat), data = forgraph, group = factor(type)) +
  geom_bar(stat = "identity", 
           position = position_dodge(width = 0.9)) +
  geom_errorbar(aes(x=type, ymin=means-error, ymax=means+error),
                width = 0.25,
                position = position_dodge(width = 0.9)) +
  xlab("") +
  ylab("") +
  scale_y_discrete(limits = c(1,2,3,4,5),
                   labels = c("Strongly\nOppose (1)", 
                              "Somewhat\nOppose (2)",
                              "Neutral (3)", 
                              "Somewhat\nFavor (4)", 
                              "Strongly\nFavor (5)"),
                   expand = expansion(add = c(0,1))) +
  scale_fill_discrete(name = "", labels = c("Denial", "Overt")) +
  theme(axis.text.y = element_text(angle = 0, size = 12.5),
        axis.text.x = element_text(size = 12.5),
        legend.text = element_text(size = 12.5)) 
ggsave("Tables&Figures\\avg_escalation_iran3.png", width = 6, height = 4, unit = "in")




